home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_4 / a4_p.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  6.3 KB  |  239 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 4.p (Pade rational Approximation).
  9. % Section    4.6,    Pade Approximations, Page 249
  10. echo on; clc; format short; hold off; clear
  11.  
  12. % Investigation of a Pade rational approximations.
  13. %                          Pn(x)
  14. %      f(x)  ≈  Rn,m(x) = -------
  15. %                          Qm(x)
  16. % where the degrees of Pn and Qm are n and m, respectively.
  17.  
  18. % The algorithm requires the Maclaurin coefficients of f(x).
  19. % Coefficient lists for several functions have been
  20. % stored in M-files named;   zcos  zsin  ztan  zexp
  21. % zacos  zasin  zatan  zcosh   zsinh   zsqrt   zlog
  22. % zsqrt4   zinv   zemx2d2   zcosde   zsinde   zlogq
  23.  
  24. % Remark. pade.m is used for Algorithm 4.p
  25.  
  26. pause    % Press any key continue.
  27.  
  28. clc;
  29. %       Pade rational approximations for  f(x) = cos(x).
  30.  
  31. % Issue the command  zcos  to load the coefficients
  32.  
  33. % into the array  C.  The function name is loaded
  34.  
  35. % into the variable fun, the degree is loaded into  m.
  36.  
  37. % The endpoints of [a,b] are loaded into  a  and b.
  38.  
  39. % Load the Taylor coefficients.
  40.  
  41. [fun,dfun,ifun,x0,m,C,Ax] = zcos;
  42.  
  43. pause % Press any key to continue.
  44.  
  45. clc;
  46.  
  47. % Place the degree of Pn(x) in  n
  48.  
  49. % Place the degree of Qm(x) in  m
  50.  
  51. % Remark. Some combinations of n and m are incompatible.
  52. %         This might result in a singular matrix
  53. %         or division by zero where Qm(x) = 0.
  54.  
  55. n = 2;
  56. m = 2;
  57. [P,Q] = pade(C,n,m);
  58.  
  59. pause % Press any key to graph f(x) and Rn,m(x). 
  60.  
  61. clc; clg;
  62. a = -3;    % You can change the left  endpoint a.
  63. b =  3;    % You can change the right endpoint b.
  64. c = -1.5;
  65. d =  1;
  66. X = a:.05:b;
  67. x = X;
  68. Y = eval(fun);
  69. axis([a b c d]);...
  70. R = polyval(P,X)./polyval(Q,X);...
  71. plot(X,Y,'-g',X,R,'-r');...
  72. hold on;...
  73. plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
  74. xlabel('x');...
  75. ylabel('y');...
  76. Mx1 = ['Comparison of ',fun,' and R'];...
  77. Mx2 = [Mx1,num2str(n),',',num2str(m),'(x)'];...
  78. title(Mx2);...
  79. grid;...
  80. axis;...
  81. hold off;...
  82. shg; pause    % Press any key to continue.
  83.  
  84. Mx1 = 'The function is  f(x) = ';
  85. Mx2 = 'The interval is  ';
  86. Mx3 = 'Pn(x) = p(1)x^n + p(2)x^(n-1) + ... + p(n)x + p(n+1)';
  87. Mx4 = 'The degree is  n = ';
  88. Mx5 = ', and the coefficients list  P  is:';
  89. Mx6 = 'Qm(x) = q(1)x^m + q(2)x^(m-1) + ... + q(m)x + q(m+1)';
  90. Mx7 = 'The degree is  m = ';
  91. Mx8 = ', and the coefficients list  Q  is:';
  92. clc,format short e,echo off,diary output,...
  93. disp(''),disp([Mx1,fun]),...
  94. disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
  95. disp(Mx3),disp([Mx4,num2str(n),Mx5]),...
  96. for i=1:5:n+1, disp(P([i:min(i+4,n+1)])); end,...
  97. disp(Mx6),disp([Mx7,num2str(m),Mx8]),...
  98. for i=1:5:m+1, disp(Q([i:min(i+4,m+1)])); end,...
  99. diary off, echo on
  100.  
  101. pause    % Press any key to graph f(x) - Rn,m(x).
  102.  
  103. clc; clg;
  104. a = -1;    % You can change the left  endpoint a.
  105. b =  1;    % You can change the right endpoint b.
  106. h = (b-a)/200;
  107. X = a:h:b;
  108. x = X;
  109. Y = eval(fun);
  110. R = polyval(P,X)./polyval(Q,X);
  111. c = min(Y-R);
  112. d = max(Y-R);
  113. axis([a b c d]);...
  114. plot(X,Y-R,'-r');...
  115. hold on;...
  116. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  117. xlabel('x');...
  118. ylabel('y');...
  119. Mx1 = ['The error;  ',fun,' - R'];...
  120. Mx2 = [Mx1,num2str(n),',',num2str(m),'(x)'];...
  121. title(Mx2);...
  122. grid;...
  123. axis;...
  124. hold off;...
  125. shg; pause    % Press any key for a list of numerical computations.
  126.  
  127. clc; format long;
  128. a = -1.2;    % You can change the left  endpoint a.
  129. b =  1.2;    % You can change the right endpoint b.
  130. X = a:0.1:b;
  131. x = X;
  132. Y = eval(fun);
  133. R = polyval(P,X)./polyval(Q,X);
  134. points = [X;Y;R;Y-R];
  135. Mx1=['Pade rational approximation of f(x) = ',fun];
  136. Mx2='     x(k)               f(x(k))            Rn,m(x(k))           error';
  137. clc,echo off,diary output,...
  138. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  139. diary off,echo on
  140.  
  141. pause % Press any key to construct a Pade rational approximation.
  142.  
  143. clc;
  144.  
  145. % Place the degree of Pn(x) in  n
  146.  
  147. % Place the degree of Qm(x) in  m
  148.  
  149. % Remark. Some combinations of n and m are incompatible.
  150. %         This might result in a singular matrix
  151. %         or division by zero where Qm(x) = 0.
  152.  
  153. n = 4;
  154. m = 4;
  155. [P,Q] = pade(C,n,m);
  156.  
  157. pause % Press any key to graph f(x) and Rn,m(x). 
  158.  
  159. clc; clg;
  160. a = -5;    % You can change the left  endpoint a.
  161. b =  5;    % You can change the right endpoint b.
  162. c = -1;
  163. d =  1;
  164. X = a:.05:b;
  165. x = X;
  166. Y = eval(fun);
  167. axis([a b c d]);...
  168. R = polyval(P,X)./polyval(Q,X);...
  169. plot(X,Y,'-g',X,R,'-r');...
  170. hold on;...
  171. plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
  172. xlabel('x');...
  173. ylabel('y');...
  174. Mx1 = ['Comparison of ',fun,' and R'];...
  175. Mx2 = [Mx1,num2str(n),',',num2str(m),'(x)'];...
  176. title(Mx2);...
  177. grid;...
  178. axis;...
  179. hold off;...
  180. shg; pause    % Press any key to continue.
  181.  
  182. Mx1 = 'The function is  f(x) = ';
  183. Mx2 = 'The interval is  ';
  184. Mx3 = 'Pn(x) = p(1)x^n + p(2)x^(n-1) + ... + p(n)x + p(n+1)';
  185. Mx4 = 'The degree is  n = ';
  186. Mx5 = ', and the coefficients list  P  is:';
  187. Mx6 = 'Qm(x) = q(1)x^m + q(2)x^(m-1) + ... + q(m)x + q(m+1)';
  188. Mx7 = 'The degree is  m = ';
  189. Mx8 = ', and the coefficients list  Q  is:';
  190. clc,format short e,echo off,diary output,...
  191. disp(''),disp([Mx1,fun]),...
  192. disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
  193. disp(Mx3),disp([Mx4,num2str(n),Mx5]),...
  194. for i=1:5:n+1, disp(P([i:min(i+4,n+1)])); end,...
  195. disp(Mx6),disp([Mx7,num2str(m),Mx8]),...
  196. for i=1:5:m+1, disp(Q([i:min(i+4,m+1)])); end,...
  197. diary off, echo on
  198.  
  199. pause    % Press any key to graph f(x) - Rn,m(x).
  200.  
  201. clc; clg;
  202. a = -1;    % You can change the left  endpoint a.
  203. b =  1;    % You can change the right endpoint b.
  204. h = (b-a)/200;
  205. X = a:h:b;
  206. x = X;
  207. Y = eval(fun);
  208. R = polyval(P,X)./polyval(Q,X);
  209. c = min(Y-R);
  210. d = max(Y-R);
  211. axis([a b c d]);...
  212. plot(X,Y-R,'-r');...
  213. hold on;...
  214. plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
  215. xlabel('x');...
  216. ylabel('y');...
  217. Mx1 = ['The error;  ',fun,' - R'];...
  218. Mx2 = [Mx1,num2str(n),',',num2str(m),'(x)'];...
  219. title(Mx2);...
  220. grid;...
  221. axis;...
  222. hold off;...
  223. shg; pause    % Press any key for a list of numerical computations..
  224.  
  225. clc; format long;
  226. a = -1.2;    % You can change the left  endpoint a.
  227. b =  1.2;    % You can change the right endpoint b.
  228. X = a:0.1:b;
  229. x = X;
  230. Y = eval(fun);
  231. R = polyval(P,X)./polyval(Q,X);
  232. points = [X;Y;R;Y-R];
  233. Mx1=['Pade rational approximation of f(x) = ',fun];
  234. Mx2='     x(k)               f(x(k))            Rn,m(x(k))           error';
  235. clc,echo off,diary output,...
  236. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  237. diary off,echo on
  238.  
  239.